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ABSTRACT 


A method of finding a first approximation to a crustal magnetiza- 
tion distribution from inversion of satellite magnetic anomaly data is 
described. Magnetization is expressed as a Fourier Series in a segment 
of spherical shell. Input to this procedure is an equivalent source 
representation of the observed anomaly field. Instability of the 
inversion occurs when high frequency noise is present in the input 
data, or when the series is carried to an excessively high wave number. 
Preliminary results are given for the United States and adjacent areas. 
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INTRODUCTION 

The polar-orbiting satellites OGO 2, 4, and 6 collected total-field 
niagnet'ic data at elevations above 400km, A preliminary anomaly data set 
was created by selecting data with minimal external field effects, and by 
subtracting a 13^^ degree spherical harmonic representation of the core 
field fit to this data subset. Regan et al (1975) published a l°-average 
representation of the anomaly field for a strip around the world between 
50°N and 50°S, and described the data reduction procedures. 

This paper is a review of a simple method for finding a first 
approximation to a crustal magnetization distribution which will pro- 
duce a field which reproduces the measured satellite field. The term 
"crust" is used loosely to mean a layer bounded by the Earth's surface 
and the Curie isotherm, and may or may not correspond to the petrologic 
crust in a given area. 

MODELING THE ANOMALY FIELD 

The anomaly data set is contaminated by noise of three main kinds: 
1) instrument noise, 2) local current effects, and 3) very long wave- 
length effects due to magnetospheric ring currents. The third effect 
has been described by Langel and Sweeney (1971). Cain and Davis (1973) 
modeled this effect as a first zonal harmonic, which they fit to 
individual satellite passes between 50°N and 50°S geomagnetic. Figure 1 
is three groups of three passes in profile form; the tracks are shown i 
Figure 2. Within each group the satelli^.a elevations are similar, and 1 
the profiles should be similar. The raw anomaly data is shown in row 1 
clearly, residual long wavelength effects are present in the individua 
profiles. Row 2 is "ring-corrected" data. The correction generally 
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improves the internal agreement, but a substantial residual remains, 
and some further correction is needed. This residual is partly respon- 
sible for the north-south elongation of anomaly contours, reflecting 
the satellite tracks, in the world map of Regan et al (1975). The result 
is effective high frequency noise in the east-west direction. In rov;s 
3 and 4 linear and quadratic functions, respectively, have been fit 
to the individual profiles and subtracted. In row 5, a first zonal 
harmonic term has been fit only over the latitude range shown. The 
internal agreement is greatly improved in each case, but seems slightly 
better for the quadratic fit; therefore, a quadratic function was fit. 
to and subtracted from each profile used in the computations described 
below. 

The data is distributed through a considerable elevation range, 
but we would like to be able to represent the field at an arbitrary 
constant elevation. For this reason, and to average out instrumental 
and transient current effects, the anomaly field was modeled by an 
equivalent source procedure. This consisted of setting out an array 
of dipoles at the Earth's surface in a 4*^ latitude-longitude grid, and 
determining a set of moments for the dipoles which would generate an 
artificial field which would make a least-squares best fit to the data; 
the mathematics is outlined in the Appendix. The dipoles were oriented 
along the direction of the main field, although to simply model the field 
this direction is not critical. The input data was limited to the 
elevation range 400-550 km. The fit of the computed field to the data 
is to a standard deviation of about 1 nT. Once the dipole moments are 
determined, the field can be computed at any elevation; a computation 
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at 450 km is shown in Figure 3, Figure 4 shows the fit of observed and 
computed fields for an arbitrary selection of profiles running between 
lO^N and 50°N in the area of Figure 3. 

The input field to the inversion procedure outlined below must be 
smooth. Since the equivalent soui^ce field fills this requirement, it, 
rather than the raw anomaly data, was so used. 

MAGNETIZATION DISTRIBUTION 

The set of magnetic moments determined in the equivalent sojrce 
computation vary irregularly, and have no particular physical signif- 
icance. One approach to developing a physically meaningful model of the 
magnetic source distribution is to seek a continuous distribution of 
magnetization in a layer of constant thickness which will give rise to a 
field which closely fits the input field. The result is a first approx- 
imation to gross magnetization variations in the magnetic crust. The 
procedure is similar to that for the equivalent source computations 
described above, but with two essential differences. First, !!*« mrces 
are 2° blocks 40 km thick, rather than dipoles. An approximate source 
function was developed for the anomSly due to such spherical prisms 
(see Appendix). Second, rather than allowing the moments of the sources 
to vary independently, their magnetizations were specified by the value 
of a double Fourier series in latitude and longitude having terms of the 
form 

A^j* (cos,sin)(2Tiix/X)' (cos,sin)(2irjy/Y). (1) 

The unknown parameters in the least-squares formulation are then the 
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constants of the series, rather than the magnetic moments of individual 
sources. Map areas 40*^ by 40° were treated individually. The equivalent 
source field, tapered to zero 4° beyond each map border, was used for 
input. The Fourier series was expressed within the extended area; thus, 
in expression (1) above X=Y=48°. 

Once the series parameters are determined, the field can be com- 
puted at arbitrary elevation. The result at 450 km is shown in Figure 5, 
which is to be compared with Figure 3. The magnetization distribution 
itself in units of emu/cc x 10^ is shown in Figure 6. 

There is a particular advantage to having the data so high above 
the sources. A particular source block has very nearly the same anomaly 
as a block twice the thickness and half the magnetization; thus, one 
can readily convert the model of magnetization variation in a layer of 
constant thickness to variations in a layer of variable thickness where 
there is Independent evidence on the thickness of the magnetic crust. 

SOURCES OF INSTABILITY 

High frequency components of the field tend to be strongly ampl- 
ified on inversion. Two examples of difficulties of this type are 
discussed below. 

The map on the left in Figure 7, a test area in the Indian Ocean, 
was made by averaging corrected data within the elevation range 400- 
550 km over 1° squares; average data generally contains high frequency 
noise contamination. The southern part of the map is in high magnetic 
latitudes, and external field noise is present along the southern 
border. The map on the right is the result of an inversion in v/hich 


these components have evidently been exaggerated, producing a charact- 
eristic cell-like structure. 
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Figure 8 shows a second kind of problem. There is an obvious 
question of how large the maximum wave number in the Fourier represen- 
tation of magnetization can be. Figure 8 shows the results of computations 
for maximum wave number 3,4,5 and 6. The maps in the top row are the 
magnetization results; shown below are the corresponding power ''?.it-.jtra. 
Note that while the low frequency part of the spectrum does not change 
very much, for m=5 and 6 increasing energy is entering the high frequency 
part. Apparently spurious effects appear in the corresponding magnetiz- 
ation maps. It would seem, then, that for an area this size the maximum 
appropriate wave number is 4, despite the fact that the fit of input to 
computed values continues to improve with expanding series, as indicated 
by the standard deviation values. The result suggests the limits of source 
resolution with this method. m=4 corresponds to a minimum wavelength of 
roughly twice the elevation of the data. The results of Figure 6 were 
computed for m=5, which corresponds to about the same minimum wavelength 
because the map area is larger. 
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APPENDIX 


Equivalent Source Field 

We seek an expression for the potimtial of a dipole located at a 
point j on the Earth's surface at some external measurement point i. 
Specify the coordinates as (rj, Oj, and (r^, 0^, <}.^), where r is 
radial distance, e is colatitude, and ^ is longitude east. 


where t.. is the distance between i and j. 

1 J 

= r? + r? -2r.r<Cosc*. , 

ij J 1 J 1 "ij ’ 

where is the central angle between i and j, and 

cosc^j = cosGj cosG^. + sinGj sino^ cos((|)^“(|)j), 


The components of nb are 

J 

(m.sini, m.cosIcosD, m.cosIsinD), 

J V J 

where I and D are inclination and declination of the dipole, taken 
to be that of the main field vector at the dipole. 


Differentiation and substitution yields 
V. .=m.{sinl{r .-r.A) + cosIcosDr.B-cosIsinDr.C)/a.? 

IJJ J» * 1 ij 

where A = cose .cose . + sin0 .sine .cos((|>.-(|. .) 

J I J I * J 

B = sine .COS0. -COSO .sine. cos(ij).-(^.) 

J 1 J 1 1 J 

C = sine^.sin(.}.^.-.j.j). 


The gradient of \t-^ in the total field direction is the anomaly in the 
total field, B... 

I J 


The anomaly due to ell the dipoles is 


Using a procedure outlined by Cain et al (1967), we determine a set 

of values for the m. which will minimize the square residuals between 

J 

observed and measured F over all points i. 


Field due to Spherical Prison of Elemental Area 

Proceeding from the dipole result, we write down an expression for 
the potential of a volume element. 


J 

where is the magnetization vector. 

r 0 <p 

Then integrating in the r direction, and replacing the infinitesimal 
angles by small finite angles gives the potential of a spherical prism 


of elemental area. 




fRz rj 
— ii., dr. 

Ri ^ij J 


•"(-Vi^ + Jgr.B.J^r.C)sin0j.a6.A.frj 


Rz 


To* dr 

Ri ^ 


where (R2-R1) = 40 km, the approximate "crustal" thickness used in this 
paper. Experiments with computations have shown that the A0 and a<{) 
can be taken to be 2° in the above expression, and still give an 
excellent approximation to the field of a prism with finite angular 
dimensions. 

The integrals above are given by Gradsbteyn and Ryzhik (1965); the 
expressions are fairly lengthy, and are not reproduced here. 
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Figure 4: Corrected measurements for arbitrarily selected group of 

profiles in the area of Figure 3 (dots) and computed values 
(pluses). Computations cease where profile passes out of 
geographic area or (more coninonly) rises above 550 km elevat'Jon. 




Figure 5: Anomaly field at 450 km computed from magnetization distributi 

of Figure 6. Units are nT. 
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Figure 8: Magnetization distributions' and corresponding power spectra 

for maximuni wave number m= 3»4,5, and 6. Scale shown in wave 
number plus one. 



